/**
  CNOK project, Anyang Normal University, IMP-CAS
  \class TARs
  \brief calculate the spectroscopic reduction factor Rs=sig_exp/sig_th for
  intermediate-energy single-nucleon knockout reactions. The main function of
  this class is to calculate the inclusive single-nucleon removal c.s.
  \author SUN Yazhou, asia.rabbit@163.com
  \since 2022/05/23
  \date Last modified: 2022/05/24 by SUN Yazhou
  \copyright 2020-2023 SUN Yazhou
  \copyright CNOK project, Anyang Normal University, IMP-CAS
*/

#include <iostream>
#include <iomanip>
#include <fstream>
#include <cmath>
#include <cstring>
#include "TARs.h"
#include "TAException.h"
#include "TAParaManager.h"

using std::cout;
using std::endl;
using std::setw;
using std::ofstream;

const map<char, int> TARs::LL = {
  {'s', 0}, {'p', 1}, {'d', 2}, {'f', 3}, {'g', 4}, {'h', 5}, {'i', 6}
};

// is_p: true: -p; false: -n
TARs::TARs(int a, double sp, double sn, double cse, double dcse) : fIsProtonRemoval(false),
    fA(a), fSp(sp), fSn(sn), fCSExp(cse, dcse){
  if(a <= 0 || a >= 400) TAException::Error("TARs", "ctor: invalid A: %d", a);
  Initialize();
} // end ctor

void TARs::YamlName(const string &orbit, char *cfg){
  int len = orbit.size();
  int pos = 3; strncpy(cfg, orbit.c_str(), 3);
  if('/' == orbit[pos]) pos += 3; // skip '/2'
  cfg[3] = orbit[pos]; pos += 1; // '[+-]' or '/'
  if('/' == orbit[pos]) pos += 2; // skip '/2'
  if('+' == orbit[pos]) cfg[4] = 'p';
  else if('-' == orbit[pos]) cfg[4] = 'm';
  else TAException::Error("TARs", "YamlName: illegal config: %s", orbit.c_str());
  pos += 2;
  // level no. for the same Jpi
  if(len > pos) cfg[5] = orbit[pos];
  strcat(cfg, ".yaml");
} // end member function YamlName

// to tell if b is a legal orbit: e.g.: 0p1/2_1/2-, or 1s1/2_1-_1, or
bool TARs::ob(const string &b){
  int p = 0; bool r = true;
  if(!isdigit(b[p++])) r = false;
  const char l[] = "spdfghi"; bool ll = false;
  for(const char c : l) if(c == b[p]) ll = true;
  if(!ll) r = false;
  if(!isdigit(b[++p])) r = false; // nominator
  if('/' == b[p+1] && '2' == b[p+2]) p += 3; // '/2'
  else r = false;
  if(b[p++] != '_') r = false;
  if(!isdigit(b[p])) r = false; // nominator
  if('/' == b[p+1] && '2' == b[p+2]) p += 2; // '/2'
  p++;
  if('+' != b[p] && '-' != b[p]) r = false; // parity
  if((++p) == int(b.size())) return r;
  if(b[p] != '_') r = false;
  return isdigit(b[++p]) && r;
} // end inline function ob

/// \param csth: theoretical c.s.
void TARs::AddConfig(const string &orbit, double ex, double c2s, double cssp){
  if(ob(orbit)) fOrbit.push_back(orbit);
  else TAException::Error("TARs", "AddConfig: illegal input of configuration: %s\
should be of the format \'nlj_Jpi[_No]\'", orbit.c_str());
  fEx.push_back(ex);
  fC2S.push_back(c2s);
  fCSsp.push_back(cssp);
} // end member function AddConfig

void TARs::Initialize(){
  fEx.clear();
  fC2S.clear();
  fCSsp.clear();
  fCSth.clear();
} // end member function Initialize

// calculate the Rs
void TARs::Rs(){
  fCSthT = 0.;
  if(fIsProtonRemoval) fDeltaS = fSp - fSn;
  else fDeltaS = fSn - fSp;
  double ExAvrg = 0.;
  const int nn = fOrbit.size();
	for(int i = 0; i < nn; i++){
		int n, l, j; char lm;
		sscanf(fOrbit[i].c_str(), "%d%c%d%*s", &n, &lm, &j);
		if(LL.count(lm)) l = LL.at(lm);
		else TAException::Error("TARs", "Rs: orbital momentum notation %c not found", lm);
		const int N = 2.*n + l;
		// const int N = 0;
		fCSth.push_back(pow(fA/(fA-1.), N)*fC2S[i]*fCSsp[i]);
		ExAvrg += fCSth[i]*fEx[i];
		fCSthT += fCSth[i];
	} // end for over i
	ExAvrg /= fCSthT;
	fDeltaS += ExAvrg;
  fRs = fCSExp / fCSthT;

	cout << std::fixed << std::setprecision(6);
	cout << "theoretical c.s.:" << endl;
	cout << std::left;
	for(int i = 0; i < nn; i++){
    cout << setw(15) << fOrbit[i] << " sigma_sp " << fCSsp[i];
    cout << " mb, sigma_th " << fCSth[i] << " mb" << endl;
  } // end for over i
	cout << "Total theoretical c.s.: " << fCSthT << " mb" << endl;
	cout << "\033[32;1mRs: " << fRs << endl;
	cout << "\033[0mDelta S: " << fDeltaS << " MeV" << endl;

  char tt[256]{}; strcpy(tt, TAException::time0(true)); // suitable for filename
  char tt1[256]; strcpy(tt1, TAException::time0()); // suitable for human reading
  char file[1024]; sprintf(file, "%ssummary_%s.txt", PARA->BaseDir(), tt);
  ofstream fout(file); fout << std::setprecision(6) << std::fixed;
  fout << "Run at " << tt1 << endl;
  fout << "theoretical c.s.:" << endl;
	fout << std::left;
	for(int i = 0; i < nn; i++){
    fout << setw(15) << fOrbit[i] << " sigma_sp " << fCSsp[i];
    fout << " mb, sigma_th " << fCSth[i] << " mb" << endl;
  } // end for over i
	fout << "Total theoretical c.s.: " << fCSthT << " mb" << endl;
	fout << "Rs: " << fRs << endl;
	fout << "Delta S: " << fDeltaS << " MeV" << endl;
  fout.close();
} // end member function Rs
